act_site = [1 2 4 7 8];
pas_site = [3 5 6 9];
con_site = 10;

for site = 1:9
    for layer = 1:4
        for time = 1:4
            cur_act = proj_meta(site).rd(layer,time).act(:,1:7500);

            % generate binary activity during feedback.
            cur_bi = proj_meta(site).rd(layer,time).velM_smoothed(1:7500);
            cur_bi = -cur_bi; % make the velocity of forward-running positive.
            cur_bi(cur_bi > 0.001) = 1;
            cur_bi(cur_bi <= 0.001) = 0;
            
            % generate the frame id of ps onset.
            ps_ons = find(diff(proj_meta(site).rd(layer,time).ps_id(1:7500))>1)+1;
            % remove early & late ps.
            ps_ons(ps_ons > 7425) = [];
            ps_ons(ps_ons < 6) = [];

            rps_nbr = 0; % number of real ps
            rps_id = [];
            for ps_id = 1:length(ps_ons)
                if sum(cur_bi(ps_ons(ps_id)-5:ps_ons(ps_id)+15)) == 21
                    rps_nbr = rps_nbr+1;
                    rps_id(rps_nbr) = ps_id;
                    all_rps_nbr(site,time) = rps_nbr;
                end
            end

            if rps_nbr >= 10
                rps_id = rps_id(1:10);
%             else
%                 rps_id = [];
            end
            
            for cell_id = 1:size(cur_act,1)
                rps_resp_cnt = 0;
                rps_resp = [];
                for rps_act_nbr = 1:length(rps_id)
                    rps_frame_id = ps_ons(rps_id(rps_act_nbr));
                    rps_resp_cnt = rps_resp_cnt+1;
                    prps_resp(rps_resp_cnt) = max(cur_act(cell_id,rps_frame_id+8:rps_frame_id+22));
                end
                prps_resp_by_cell{site,layer,time,cell_id} = prps_resp;
            end
        end
        cell_cnt{site,layer} = size(proj_meta(site).rd(layer,1).act,1);
    end
end

for site = 1:9
    for layer = 1:4
        if cell_cnt{site,layer} > 0
            for cell_id = 1:cell_cnt{site,layer}
                cur_cell_act = horzcat(prps_resp_by_cell{site,layer,:,cell_id});
                mprps{site,layer,cell_id} = mean(cur_cell_act);
            end
        end
    end
end